Collective modes of coupled phase oscillators with delayed coupling 
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We study the effects of delayed coupling on timing and pattern formation in spatially extended 
systems of dynamic oscillators. Starting from a discrete lattice of coupled oscillators, we derive a 
generic continuum theory for collective modes of long wavelength. We use this approach to study 
spatial phase profiles of cellular oscillators in the segmentation clock, a dynamic patterning system 
of vertebrate embryos. Collective wave patterns result from the interplay of coupling delays and 
moving boundary conditions. We show that the phase profiles of collective modes depend on coupling 
delays. 
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In complex dynamical systems interactions between dif- 
ferent elements can give rise to dynamical order and 
spatio-temporal patterning ■ These interactions can 
themselves be the result of a complex process. Therefore 
coupling can have internal dynamics that involve time 
delays. The importance of time delays in the coupling 
of oscillators was first recognized by Schuster and Wag- 
ner Q , who showed that two oscillators can entrain even 
if coupling is delayed. In this case, multistability of dy- 
namic states can occur as a consequence of time delays 
and the collective frequency of the system depends on 
the delay time. The roles of time delays have been ad- 
dressed in studies of different oscillator systems. It was 
found that coupling delays can give rise to a rich variety 
of behaviors [HPTl^] . 

Significant time delays in the coupling of oscillators oc- 
cur in many systems in biology, engineering and physics. 
Coupling delays between a sending and a receiving ele- 
ment can arise due to (i) intrinsic times of signal gener- 
ation in the sending element, (ii) the finite propagation 
velocity of signals, and (iii) the slow signal processing of 
the receiving clement. Coupling delays of type (ii) are in- 
evitable in some engineered systems [13l4l7l | and neuronal 
systems [l8- 21|. In some cases such delays can be used 
to implement control schemes [13, H3| ■ In the context of 
signaling processes between biological cells 2J-|26(, cou- 
pling delays of type (i) and (iii) occur naturally because 
of the complex internal kinetics of intercellular signaling. 

An important biological example in which the delayed 
coupling of dynamic oscillators plays a key role for the 
formation of patterns is the the so-called segmentation 
clock [27], [28| . This system operates during embryonic 
development of all vertebrate animals, generating a seg- 
mented morphology along the vertebrate body axis of the 
embryo. These segments, called somites, are the embry- 
onic precursors of adult vertebrae. Somites are formed 
sequentially in a dynamic tissue which elongates during 
the process. The segmentation clock is thus a rhyth- 



mic pattern generator in the tissue resulting from the 
collective organization of many cells. Each cell is an au- 
tonomous oscillator with time-periodic activation of cer- 
tain genes 0, H^] • These oscillators are noisy and are 
coordinated by intercellular signaling (30j . 

It was recently suggested that coupling delays between 
genetic oscillators play an important role for the dynam- 
ics of the segmentation process and that they influence 
the collective frequency of cellular oscillations [3lJ. By 
comparing theory with quantitative experiments in fish 
embryos it was subsequently shown that coupling delays 
indeed play a crucial role in pattern formation and that 
the collective oscillation frequency is altered in fish with 
mutations affecting the coupling process between oscil- 
lators 32|. The segmentation clock is therefore a prime 
example of a population of coupled oscillators in which 
coupling delays play a crucial role. Its understanding re- 
quires a theory of locally coupled oscillators with delays 
in a spatially extended system. 

Despite their general interest, the effects of coupling 
delays in spatially extended systems are still poorly un- 
derstood. It has been shown, using a low dimensional 
approximation [33| , that such systems can display a wide 
range of spatiotemporal patterns [34[ . An important ap- 
proach to describe spatially extended systems is to fo- 
cus on long-wavelength modes, which can be described 
by a simplified continuum limit which ignores details on 
small scales. Continuum descriptions have been devel- 
oped for chains of coupled phase oscillators without delay 
Mil and for rings of coupled phase oscillators in the 
limit of short delay times [lfj . It was shown that short 
time delays are equivalent to a phase shift in the coupling 
[37l . |38| . This simplification is lost when coupling delays 
are longer than the time scale defined by the coupling 
strength. In this case time delays have to be considered 
explicitly [H]. 

In this Letter, we study the collective oscillatory modes 
of the vertebrate segmentation clock. We introduce a 
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generic continuum description for the long wavelength 
modes of spatially extended lattices of coupled oscillators 
with arbitrary coupling delays. We use this theory to 
discuss the emergent collective frequency, as well as phase 
profiles of collective modes. 

The vertebrate segmentation clock is a tissue level pat- 
tern generator [27|. It consists of a population of cells 



that collectively generate oscillating activity patterns of 
cyclic genes in the tissue. Neighboring cells are cou- 
pled via a slow molecular signaling system that intro- 
duces time delays [31, 32 1. Thus, we consider a system 
of coupled oscillators defined on a regular d-dimensional 
hypercubic lattice with lattice spacing a. Oscillators sit 
on lattice sites with position Xj = ai with i = (ii, ..,id), 
i n G Z. The evolution equation for the phase 0i of the 
oscillator i, coupled to its nearest neighbors j, is 



d9j{t) 
dt 
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hm - r) - 9i(t)) , (1) 



where u>i is the intrinsic frequency of the oscillator i, £j de- 
notes the coupling strength of this oscillator to its neigh- 
bors, and r > is a time delay in the coupling. The 
coupling is described by the 27r-periodic function h. 

We are interested in long wavelength collective modes, 
for which the phases 9\ vary smoothly over distances that 
are long compared to the lattice spacing a 35|. We as- 



sume that intrinsic frequency and coupling vary smoothly 
in space so as to fulfill this condition. In this situation, 
we can approximate the phases of the oscillators in the 
discrete system by a continuum phase field 0(x, i). In 
the following we derive an equation for 0(x, t), which 
describes the long wavelength dynamics of the discrete 
oscillator system, Eq. (JXJ) . To simplify the notation, we 
derive the continuum theory on a one-dimensional lat- 
tice. The generalization to higher dimensions is straight- 
forward. We first perform a Taylor expansion of the cou- 
pling function h in Eq. ([I]) in powers of a, assuming 
that a/t <C 1, where I ~ (dO/dx)" 1 is a characteristic 
wavelength of the collective mode, 
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The prime denotes the derivative of h with respect to its 
argument, and we have defined 

A±(x,t)=6(x±a,t-T)-6(x,t), (3) 
A T (x,t) = 6(x,t-r) -0(x,t), (4) 
6 T {x,t) =9(x,t-r). (5) 

Introducing continuum functions tj(x) and e(x) of the in- 
trinsic frequency and coupling strength of the oscillators, 
together with Eq. @, we obtain a continuum description. 



The phase field obeys 
d6(x,t) 



dt 



■ u(x) + e(x)h(A T (x,t)) 
+ £ -^h"(A T (x,t)) 



dO{x,t-r) 



(0) 



e(x)a 2 



d 2 0(xA 



+ ^-h'(A T (x,t))- 

Note that the first order contribution in a vanishes. 
Higher order terms can be neglected in the long wave- 
length limit and have been dropped here. 
For arbitrary dimension we find 



d0(x,t) 



dt 



= w(x) + e{x)h(9(x, t-r)- 0(x, t)) 
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h" (0(x, t-r) - 0(x, t)) [V0(x, t - t)] 2 
ti (0(x, t-r) - 0(x, t)) V 2 0(x, t-r) , 
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where x G H d is a position vector in d-dimensional space. 
The coupling delay r enters in the spatial derivatives of 
the phase field, as well as in the arguments of the cou- 
pling function and its derivatives. For a vanishing delay, 
t = 0, we recover the classical case of locally coupled 
oscillators (35|. Eq. (J7J describes the collective modes 
of general extended systems of oscillators with coupling 
delays. 

We now employ our approach to describe gene activity 
patterns of cells in the segmentation clock of vertebrate 
embryos. We simplify our description of the segmen- 
tation process by using a semi-infinite one-dimensional 
geometry, Fig. [U in which a system of oscillators is de- 
scribed by Eq. ©. The system in which oscillators cre- 
ate patterns extends from the anterior for x — > — oo to 
a posterior boundary at x p (t) — xq + vt. This moving 
boundary describes the elongating tip of the tissue. Here 
v is the elongation velocity of the tissue and v/a is the 
rate at which cells are added at the extending end. In 
the following we choose xo = without loss of generality. 

In the segmentation clock, concentration gradients of 
signaling molecules exist across the tissue. The source 
of these gradients is located at the tip of the elongating 
tissue, and the resulting concentration gradients move 
together with the tip. Perturbations of these molecu- 
lar gradients produce effects consistent with alteration of 
the intrinsic frequency of the cellular oscillators 3^, 40j |. 
Motivated by these observations, the spatial profiles uj(x) 
and e{x) are assumed to depend only on the distance d = 
x p — x from the elongating tip [3ll ItH - II^ , i.e. they travel 
together with the moving boundary. Therefore, we intro- 
duce a reference frame co-moving with the tip boundary, 
where the frequency and coupling strength profiles are 
stationary. We define the coordinate y = x — vt and 



the co-moving phase field $(y, t) = 6{y + vt,t). Since 
d = x p — x — —y, the functions e(y) and ui(y) are time- 
independent. We assume that these functions take max- 
imal values w(0) = ui and e(0) = £o at y = 0, that they 
decay monotonically and vanish in the limit y — > — oo. 

The fact that both e(y) and ui(y) vanish for large neg- 
ative y implies, according to Eq. (|6]), that the pattern 
0(x, t) becomes stationary for x — > — oo, i.e. dO/dt ~ in 
this limit. This stationary pattern describes the develop- 
ing segments. Pattern forming solutions can be obtained 
by the ansatz: 



frequency profile 



homogeneous 
frequency 



i%,t) = 0(y)+fit. 



(8) 



Using this expression, the differential equation for the 
phase profile reads for y < 

ft = v<f/(y) + uj(y) + e{y)h{<P{y + vr) - <%) - Or) 

+ £i^V {4>{y + vt) _ ^y) _ n T ) {<j)'{y + vr)) 2 



+ E^K h ' + VT ) _ 0( y ) _ Qr) <j>"{y + vr) 



(9) 



The first term in Eq. © is a drift describing phase trans- 
port due to the motion of oscillators in the co-moving 
system. The moving boundary and the co-moving pro- 
files of e and u, together with coupling delays, introduce 
non-local effects in Eq. (|9|). 

Eq. ((9]) is solved imposing boundary conditions at the 
moving tip, y = 0. Because of the nonlocal terms in- 
volving oscillators at y + vt in Eq. ([9]), the boundary 
values of 0(0), (f>'(0) and <fi"(0) arc not sufficient to de- 
termine the solution. Rather, it is necessary to specify 
the function (j)(y) in the interval < y < vt. This nonlo- 
cal boundary condition reflects the fact that the history 
of oscillators that are attached to the end is important. 
This history is specified by this boundary condition. We 
choose 4>(y) = <f>o for y > as well as <f>'(0) — and 
4>"{0) = 0. This choice corresponds to the assumption 
that all oscillators that enter the system at the tip os- 
cillate with intrinsic frequency ujq and are coupled with 
strength Eq, being in phase with their neighbors, Fig. [1] 
(red). 

Eq. §9§ together with these boundary conditions imply 
that the collective frequency ft obeys 



Q = ojq + Soh (— fir) 



(10) 



This transcendental equation is known to determine the 
collective frequency Q in a population of identical oscil- 
lators. It can have one or several coexisting solutions for 
the collective frequency 0, @, |45[ . 

The solutions of Eq. (|9]) with the boundary conditions 
imposed here, describe a dynamic biological pattern gen- 
erator that oscillates with collective frequency dd / 'dt = Q 
and produces a stationary spatially periodic structure 




OOO0000000 

Xp(t) = V t 
posterior boundary 



anterior 



FIG. 1: One dimensional lattice of coupled oscillators repre- 
senting the segmentation clock with a frequency profile (blue) 
and posterior boundary at position x v moving with velocity 
v (red). The oscillators self-organize in a phase pattern. The 
frequency profile starts at the posterior boundary with value 
u)q and drops to zero as x — > — oo. 



of length S described by sin(6'(a;)) ~ sin(27ra;/S'). It is 
formed by the collective mode of coupled oscillators that 
corresponds to traveling waves induced at the moving 
tip that propagate towards the anterior where they slow 
down and stop. From Eq. ([5]), we find 4>{y) ~ yQ/v in the 
limit y — > — oo in which both u> and e vanish, implying 



S = vT , 



(11) 



where T = 2n/Q, is the collective period. This result 
shows that we recover the basic property of a general 
clock and wavefront mechanism of vertebrate segmenta- 
tion [27, 46 1 . The segment length S depends only on the 
elongation velocity and the collective frequency of the 
oscillators at the posterior boundary. It is independent 
of the shapes of the profiles of frequency and coupling 
strength. 

Solutions to Eq. (|9|) for boundary conditions given 
above are displayed in Fig. [2] (lines) for two different 
time delays. As an example, realistic parameters pre- 
viously determined for the zebrafish segmentation clock 
have been usedjjncluding typical profiles of the functions 
e(x) and uj(x) [32j. These profiles and parameter values 
are given in the caption of Fig. [2] These solutions are 
compared to simulations of the discrete oscillator sys- 
tem, Eq. ([U, with the same parameters (dots). The 
comparison shows that our continuum theory provides 
an excellent approximation for the behavior of the dis- 
crete system. We have checked that small perturbations 
to the discrete patterns relax back to the time-periodic 
states shown in Fig. indicating that these patterns are 
locally stable. 

Fig. [5] also highlights the role of coupling delays in 
shaping the phase patterns. The solid black line shows 
the phase profile for a time delay of r = 20.75 min [32| . 
while the dashed line was obtained for r = 44.25 min 
keeping other parameters the same. Thus, the value of 
r in both cases differs by the collective period T = 23.5 
min. Note that Eq. (|10[) . which determines T, is invariant 
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our theory predicts through Eq. (fT0[) the dependence of 
the collective frequency on coupling str eng th and cou- 
pling delay, consistent with experiments [32 1 . 
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FIG. 2: Phase profiles of coupled oscillators with delayed cou- 
pling for the segmentation oscillator. Upper panel: plots of 
sin cj>(y) corresponding to the phase profiles in the lower panel, 
illustrating cyclic gene expression patterns in vertebrate seg- 
mentation. Both simulations have the same collective period 
T and elongation velocity v, making the resulting segment 
length S equal in both. Lower panel: steady state phase pro- 
files 4>{y) in the co-moving frame obtained from the continuum 
approximation Eq. (lines), and numerical simulations of 
the discrete model given by Eq. |T} (dots). We use parameters 
previously determined for the zebrafish segmentation clock 
[33 ] . resulting in a collective period T = 23.5 min for both 
cases. The coupling function is h(A) = sin(A). Boundary 
condition is 4>{y)\v>o = const. The time units are minutes 
and the unit of length a corresponds to one cell diameter, 
a ~ 10 j-im. Only every second discrete oscillator is dis- 
played, and the arrested oscillators to the left of y = —39 are 
not shown. Parameters are: v — 0.249 cell diameters/min, 
e(y) — So for y > —39, with eo = 0.07 min -1 , e(y) = min - 
for y < -39, u(y) = wo(l - e -<"+ 39 >/ 27 )/(l - e -39 / 27 ) for 
-39 < y < 0, with u = 0.2205 min -1 , u\y) = min -1 for 
y < —39 and to(y) — luo for y > 0. The biolo gica l motivation 
of the choice of e(y) and u(y) is discussed in [311 ]. 



under the transformation r — > r + mT for any integer 
to. As a consequence, the two systems shown in Fig. 
[5] oscillate with the same period. However, the phase 
profiles of these solutions differ because Eq. §§§ is not 
invariant under this transformation. This effect may offer 
a way to distinguish between different coupling delays 
producing the same period in the segmentation clock, by 
means of phase profile data. The shape of the observed 
pattern does not change when we introduce a white noise 
term to the discrete model, Eq. (JTJ) . 

Our theory can account for several observed proper- 
ties of the segmentation clock. First, we showed that the 
simple clock and wavefront relation S = vT 2?], 4(| 47 1 
between segment length S, elongation velocity v and os- 
cillation period T, holds under very general conditions of 
frequency profiles, time delays and coupling functions. 
Second, our theory produces phase patterns that can 
quantitatively fit gene expression patterns observed in 
the segmentation clock [3l|, [32|, see Fig. [5J Third, in- 
dependently of the details of the frequency and coupling 
strength profiles and the shape of the coupling function, 



Delays in oscillator coupling play an important role for 
the stability of synchronized phase profiles, as has been 
shown for simple cases [IB]]. Stable solutions of Eq. ([7]) 
describe long wavelength modes of the coupled oscillator 
system given by Eq. (JTJ) . The stability of these modes 
depends on time delays and boundary conditions. Insta- 
bilities of Eq. ([7]) respect to short wavelengths describe 
situations where the continuum description breaks down. 

Systems of coupled oscillators can display individual 
variability (iil . |49| and be subject to dynamic fluctua- 
tions [50J. Individual variability can be described by 
quenched disorder in the parameters, such as the fre- 
quency or coupling [48| . Dynamic fluctuations can be de- 
scribed by an additional noise term in Eq. ([7]) [Hl[ . With 
the addition of this noise term, Eq. ([TJ is a generaliza- 
tion of the Kardar-Parisi-Zhang equation |52j to time- 
delayed and inhomogeneous systems. In the Kardar- 
Parisi-Zhang equation the interplay of noise and non- 
linearity induces a rich phenomenology. This indicates 
that the addition of noise to our problem may have in- 
teresting effects on the dynamics. Finally, fluctuations in 
the coupling delays could be straightforwardly accounted 
for in the theory, extending Eq. ([7]) to include distributed 
coupling delays 53 1. 



In this Letter, we have introduced a continuum descrip- 
tion of long wavelength modes in extended systems of os- 
cillators with coupling delays, Eq. ([7]). We have applied 
this continuum description to a problem from biology, 
the segmentation clock. This problem involves moving 
boundaries, which together with time delays give rise to 
nonlocal effects, see Eq. ([9]) . We have proposed here that 
time delays have a key role in shaping the pattern of spa- 
tial phase profiles, apart from setting the period of the 
segmentation clock. The biological example of the seg- 
mentation clock shows that the continuum description is 
a powerful method to study extended systems with cou- 
pling delays. 
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